Contributed by: Benoît Legat
using DynamicPolynomials
using SumOfSquares
import CSDP
Consider the problem of finding both the minimum value of p = x^4 - 4x^3 - 2x^2 + 12x + 3
as well as its minimizers.
We can use SumOfSquares.jl to find such these values as follows. We first define the polynomial using DynamicPolynomials.
@polyvar x
p = x^4 - 4x^3 - 2x^2 + 12x + 3
Secondly, we create a Sum-of-Squares program searching for the maximal lower bound σ
of the polynomial.
model = SOSModel(CSDP.Optimizer)
@variable(model, σ)
@constraint(model, cref, p >= σ)
@objective(model, Max, σ)
Thirdly, solve the program and find σ = -6
as lower bound:
optimize!(model)
solution_summary(model)
Success: SDP solved Primal objective value: 0.0000000e+00 Dual objective value: 0.0000000e+00 Relative primal infeasibility: 3.29e-17 Relative dual infeasibility: 5.00e-11 Real Relative Gap: 0.00e+00 XZ Relative Gap: 1.06e-10 DIMACS error measures: 4.93e-17 0.00e+00 5.00e-11 0.00e+00 0.00e+00 1.06e-10 CSDP 6.2.0 Iter: 0 Ap: 0.00e+00 Pobj: 0.0000000e+00 Ad: 0.00e+00 Dobj: 0.0000000e+00 Iter: 1 Ap: 8.10e-01 Pobj: -1.3830330e+01 Ad: 7.29e-01 Dobj: -8.2765454e+00 Iter: 2 Ap: 9.00e-01 Pobj: -1.5490537e+01 Ad: 9.00e-01 Dobj: -1.3831579e-01 Iter: 3 Ap: 1.00e+00 Pobj: -7.5281681e+00 Ad: 9.00e-01 Dobj: -3.8241528e+00 Iter: 4 Ap: 1.00e+00 Pobj: -6.0076491e+00 Ad: 9.00e-01 Dobj: -5.7732642e+00 Iter: 5 Ap: 1.00e+00 Pobj: -6.0000210e+00 Ad: 1.00e+00 Dobj: -5.9999957e+00 Iter: 6 Ap: 1.00e+00 Pobj: -6.0000045e+00 Ad: 1.00e+00 Dobj: -6.0000015e+00 Iter: 7 Ap: 1.00e+00 Pobj: -6.0000003e+00 Ad: 1.00e+00 Dobj: -6.0000001e+00
* Solver : CSDP * Status Result count : 1 Termination status : OPTIMAL Message from the solver: "Problem solved to optimality." * Candidate solution (result #1) Primal status : FEASIBLE_POINT Dual status : FEASIBLE_POINT Objective value : -6.00000e+00 Dual objective value : -6.00000e+00 * Work counters Solve time (sec) : 1.00708e-03
We can look at the certificate that σ = -6
is a lower bound:
sos_dec = sos_decomposition(cref, 1e-4)
(-3.000000004964362 - 1.999999998797403*x + 0.9999999995176428*x^2)^2
Indeed, p + 6 = (x^2 - 2x - 3)^2
so p ≥ -6
.
We can now find the minimizers from the moment matrix:
ν = moment_matrix(cref)
ν.Q
3×3 SymMatrix{Float64}: 1.0 0.0669168 3.13383 0.0669168 3.13383 6.46842 3.13383 6.46842 22.3383
This matrix is the convex combination of the moment matrices corresponding to two atomic measures at -1
and 3
which allows us to conclude that -1
and 3
are global minimizers.
η = atomic_measure(ν, 1e-4)
minimizers = [η.atoms[1].center; η.atoms[2].center]
2-element Vector{Float64}: -1.0000001975826238 2.9999999542477047
Below are more details on what we mean by convex combination. The moment matrix of the atomic measure at the first minimizer is:
η1 = moment_matrix(dirac(monomials(x, 0:4), x => round(minimizers[1])), ν.basis.monomials)
η1.Q
3×3 SymMatrix{Float64}: 1.0 -1.0 1.0 -1.0 1.0 -1.0 1.0 -1.0 1.0
The moment matrix of the atomic measure at the second minimizer is:
η2 = moment_matrix(dirac(monomials(x, 0:4), x => round(minimizers[2])), ν.basis.monomials)
η2.Q
3×3 SymMatrix{Float64}: 1.0 3.0 9.0 3.0 9.0 27.0 9.0 27.0 81.0
And the moment matrix is the convex combination of both:
Q12 = η1.Q * η.atoms[1].weight + η2.Q * η.atoms[2].weight
3×3 Matrix{Float64}: 1.0 0.0669169 3.13383 0.0669169 3.13383 6.46842 3.13383 6.46842 22.3383
Another way to see this (by linearity of the expectation) is that ν
is the moment matrix
of the convex combination of the two atomic measures.
The monomial basis used by default can leave a problem quite ill-conditioned for the solver. Let's try to use another basis instead:
model = SOSModel(CSDP.Optimizer)
@variable(model, σ)
@constraint(model, cheby_cref, p >= σ, basis = ChebyshevBasisFirstKind)
@objective(model, Max, σ)
optimize!(model)
solution_summary(model)
Iter: 8 Ap: 9.00e-01 Pobj: -6.0000000e+00 Ad: 1.00e+00 Dobj: -6.0000000e+00 Success: SDP solved Primal objective value: -6.0000000e+00 Dual objective value: -6.0000000e+00 Relative primal infeasibility: 6.64e-12 Relative dual infeasibility: 3.14e-10 Real Relative Gap: 2.23e-09 XZ Relative Gap: 3.25e-09 DIMACS error measures: 7.25e-12 0.00e+00 5.34e-10 0.00e+00 2.23e-09 3.25e-09 CSDP 6.2.0 Iter: 0 Ap: 0.00e+00 Pobj: 0.0000000e+00 Ad: 0.00e+00 Dobj: 0.0000000e+00 Iter: 1 Ap: 9.00e-01 Pobj: -1.8383751e+01 Ad: 7.29e-01 Dobj: -1.0343874e+01 Iter: 2 Ap: 9.00e-01 Pobj: -1.5927039e+01 Ad: 9.00e-01 Dobj: -2.3550363e+00 Iter: 3 Ap: 9.00e-01 Pobj: -7.7938220e+00 Ad: 9.00e-01 Dobj: -3.9182690e+00 Iter: 4 Ap: 9.00e-01 Pobj: -6.1761570e+00 Ad: 9.00e-01 Dobj: -5.7791841e+00 Iter: 5 Ap: 9.00e-01 Pobj: -6.0169870e+00 Ad: 9.00e-01 Dobj: -5.9771917e+00 Iter: 6 Ap: 9.00e-01 Pobj: -6.0016653e+00 Ad: 1.00e+00 Dobj: -5.9999542e+00 Iter: 7 Ap: 1.00e+00 Pobj: -6.0002842e+00 Ad: 1.00e+00 Dobj: -5.9997933e+00 Iter: 8 Ap: 9.00e-01 Pobj: -6.0000531e+00 Ad: 1.00e+00 Dobj: -5.9999938e+00 Iter: 9 Ap: 1.00e+00 Pobj: -6.0000089e+00 Ad: 1.00e+00 Dobj: -5.9999987e+00 Iter: 10 Ap: 1.00e+00 Pobj: -6.0000012e+00 Ad: 1.00e+00 Dobj: -6.0000000e+00 Iter: 11 Ap: 9.00e-01 Pobj: -6.0000002e+00 Ad: 1.00e+00 Dobj: -6.0000000e+00
* Solver : CSDP * Status Result count : 1 Termination status : OPTIMAL Message from the solver: "Problem solved to optimality." * Candidate solution (result #1) Primal status : FEASIBLE_POINT Dual status : FEASIBLE_POINT Objective value : -6.00000e+00 Dual objective value : -6.00000e+00 * Work counters Solve time (sec) : 1.46699e-03
Although the gram matrix in the monomial basis:
g = gram_matrix(cref)
@show g.basis
g.Q
g.basis = MonomialBasis([1, x, x²])
3×3 SymMatrix{Float64}: 9.0 6.0 -3.0 6.0 4.0 -2.0 -3.0 -2.0 1.0
looks different from the gram matrix in the Chebyshev basis:
cheby_g = gram_matrix(cheby_cref)
@show cheby_g.basis
cheby_g.Q
cheby_g.basis = ChebyshevBasisFirstKind([1.0, x, -1.0 + 2.0x²])
3×3 SymMatrix{Float64}: 6.25 5.0 -1.25 5.0 4.0 -1.0 -1.25 -1.0 0.25
they both yields the same Sum-of-Squares decomposition:
cheby_sos_dec = sos_decomposition(cheby_cref, 1e-4)
(-3.000000002974674 - 1.9999999995036628*x + 0.9999999997868867*x^2)^2
The gram matrix in the Chebyshev basis can be understood as follows. To express the polynomial $-x^2 + 2x + 3$ in the Chebyshev basis, we start by substituting $x$ into $\cos(\theta)$ to obtain $-\cos(\theta)^2 + 2\cos(\theta) + 3$. We now express it as a combination of $\cos(n\theta)$ for $n = 0, 1, 2$: $-(2\cos(\theta) - 1) /2 + 2 \cos(\theta) + 5/2.$ Therefore, the coefficients in the Chebyshev basis is:
cheby_coefs = [5/2, 2, -1/2]
3-element Vector{Float64}: 2.5 2.0 -0.5
We can indeed observe that we obtain the same matrix as cheby_g.Q
cheby_coefs * cheby_coefs'
3×3 Matrix{Float64}: 6.25 5.0 -1.25 5.0 4.0 -1.0 -1.25 -1.0 0.25
This notebook was generated using Literate.jl.